function [Tp,Rp,Ts,Rs]=f_2mlyr(n0,n1,k1,n2,k2,n3,d1,d2,m,fi0,lambda)

%calculate transmittance and reflectance of a multi-layer system

%Written by Feiran Huang 05/31/99. Last Update 08/13/2002 Amplitudes of the fields added.

%-------------------------

% 





for n=1:length(lambda)

   fi1=asin(sin(fi0)*n0(n)/(n1(n)-i*k1(n)));
phi1(n)=fi1;
   fi2=asin(sin(fi0)*n0(n)/(n2(n)-i*k2(n)));
phi2(n)=fi2;
   fi3=asin(sin(fi0)*n0(n)/n3(n));
phi3(n)=fi3;


   [r01p,t01p,r01s,t01s]=f_rtamp(n0(n),0,n1(n),k1(n),fi0,fi1);
   [r12p,t12p,r12s,t12s]=f_rtamp(n1(n),k1(n),n2(n),k2(n),fi1,fi2);
   [r21p,t21p,r21s,t21s]=f_rtamp(n2(n),k2(n),n1(n),k1(n),fi2,fi1);
   [r23p,t23p,r23s,t23s]=f_rtamp(n2(n),k2(n),n3(n),0,fi2,fi3);

   delta01=0;

   delta12=2*pi*(n1(n)-i*k1(n))*d1*cos(fi1)/lambda(n);
   delta21=2*pi*(n2(n)-i*k2(n))*d2*cos(fi2)/lambda(n);
   delta23=2*pi*(n2(n)-i*k2(n))*d2*cos(fi2)/lambda(n);

   C01p=[1 r01p; r01p 1];

   C12p=[exp(i*delta12) r12p*exp(i*delta12); r12p*exp(-i*delta12) exp(-i*delta12)];

   C21p=[exp(i*delta21) r21p*exp(i*delta21); r21p*exp(-i*delta21) exp(-i*delta21)];

   C23p=[exp(i*delta23) r23p*exp(i*delta23); r23p*exp(-i*delta23) exp(-i*delta23)];

   

   Ap=C01p*C12p;

   tp=t01p*t12p;

   for j=1:m-1

      Ap=Ap*C21p*C12p;

      tp=tp*t21p*t12p;

   end

   Ap=Ap*C23p;

   tp=tp*t23p;

   Rp(n)=abs(Ap(2,1)/Ap(1,1))^2;

   Tp(n)=n3(n)*cos(fi3)/n0(n)/cos(fi0)*abs(tp/Ap(1,1))^2;
   
   
   C01s=[1 r01s; r01s 1];

   C12s=[exp(i*delta12) r12s*exp(i*delta12); r12s*exp(-i*delta12) exp(-i*delta12)];

   C21s=[exp(i*delta21) r21s*exp(i*delta21); r21s*exp(-i*delta21) exp(-i*delta21)];

   C23s=[exp(i*delta23) r23s*exp(i*delta23); r23s*exp(-i*delta23) exp(-i*delta23)];

   

   As=C01s*C12s;

   ts=t01s*t12s;

   for j=1:m-1

      As=As*C21s*C12s;

      ts=ts*t21s*t12s;

   end

   As=As*C23s;

   ts=ts*t23s;

   Rs(n)=abs(As(2,1)/As(1,1))^2;

   Ts(n)=n3(n)*cos(fi3)/n0(n)/cos(fi0)*abs(ts/As(1,1))^2;
   

end












